#clean up
rm(list=ls())
options(digits=22)

#power plant documentation (eGRID dataset)
#http://www.epa.gov/cleanenergy/energy-resources/egrid/index.html

#set working directory
setwd('/Volumes/MONOGAN/polluterPlacement/data/')
#setwd('/Users/jamie/Desktop/replication/')

#load libraries
library(spatstat)
library(maptools)
library(rgdal)
library(maps)
library(foreign)
library(mgcv)
#library(rgeos)
library(spdep)
library(ggmap)
library(car)
library(CircStats)
#install.packages("/Volumes/MONOGAN/CircSpatialV3.tar.gz",repos=NULL, type="source")
library(CircSpatial)
source("altKrigCRF.R")
checkCRSArgs <- function(uprojargs) {.Call("checkCRSArgs", uprojargs, PACKAGE="rgdal")}

#options
options(scipen=8)

#Necessary function: define the haversine formula (Banerjee, Carlin & Gelfand 2004, 17-19)
haversine<-function(lon1,lon2,lat1,lat2){
	R<-6371
	x1<-(lon1*pi)/180
	x2<-(lon2*pi)/180
	y1<-(lat1*pi)/180
	y2<-(lat2*pi)/180
	dist<-R*acos((sin(y1)*sin(y2))+(cos(y1)*cos(y2)*cos(x1-x2)))
	return(dist)
	}
	
	
############STEP 1: MERGE POWER PLANT DATA WITH CONTROL GROUPS###################
#load power plant data
data.0<-read.dta('egrid_allplants.dta') # plants

#subset power plant data
data.2<-subset(data.0,subset=pstatabb!="DC"&pstatabb!="AK"&pstatabb!="HI",select=c(pname,lat,lon,pstatabb,orispl,online_year)) #Lower 48 only.
data.2$air<-1
data.2$tsdf<-0
data.2$lqg<-0
dimnames(data.2)[[2]]<-c('fac_name','latitude','longitude','loc_state','reg_id','online_year','air','tsdf','lqg')

#load and subset waste treatment site control data (cleaned already)
control.0<-read.csv('majorAirAllDist.csv')
control.1<-subset(control.0,subset=air==0)
control.1$online_year<-NA 
control.2<-subset(control.1,select=c(fac_name,latitude,longitude,loc_state,reg_id,online_year,air,tsdf,lqg))
#rm(data.0,control.0)

#append data sets together
full.0<-rbind(control.2,data.2)

############STEP 2: DATA CLEANING###################
#eliminate miscodes for latitude and longitude not being in the same state
#code for finding out-of-bounds observations
#print(attr(SD.ppp,"rejects")$x,digits=22)
#print(attr(SD.ppp,"rejects")$y,digits=22)
#> print(attr(SD.ppp,"rejects")$x,digits=22)
#[1] -96.21970367431640625
#> print(attr(SD.ppp,"rejects")$y,digits=22)
#[1] 45.221401214599609375
#full.0[full.0$loc_state=="SD",2:3]

full.0$reject<-1
full.0$reject[full.0$longitude==-85.152801513671875 & full.0$latitude==32.80030059814453125 & full.0$loc_state=="AL"]<-0
full.0$reject[full.0$longitude==-85.09110260009765625 & full.0$latitude==32.6631011962890625 & full.0$loc_state=="AL"]<-0
full.0$reject[full.0$longitude==-121.1750030517578125 & full.0$latitude==44.100299835205078125 & full.0$loc_state=="CA"]<-0
full.0$reject[full.0$longitude==-114.7003021240234375 & full.0$latitude==32.73889923095703125 & full.0$loc_state=="CA"]<-0
#full.0$reject[full.0$longitude==-86.40059661865234375 & full.0$latitude==36.315601348876953125 & full.0$loc_state=="TN"]<-0
full.0$reject[full.0$longitude==-117.11669921875 & full.0$latitude==30.666700363159179688 & full.0$loc_state=="CA"]<-0
full.0$reject[full.0$longitude==-118.153900146484375 & full.0$latitude==32.043300628662109375 & full.0$loc_state=="CA"]<-0
full.0$reject[full.0$longitude==-118.221099853515625 & full.0$latitude==33.32080078125 & full.0$loc_state=="CA"]<-0
full.0$reject[full.0$longitude==-107.96669769287109375 & full.0$latitude==36.799999237060546875 & full.0$loc_state=="CO"]<-0
full.0$reject[full.0$longitude==-90.06670379638671875 & full.0$latitude==42.06890106201171875 & full.0$loc_state=="IA"]<-0
full.0$reject[full.0$longitude==-90.74500274658203125 & full.0$latitude==40.4906005859375 & full.0$loc_state=="IA"]<-0
full.0$reject[full.0$longitude==-91.3719024658203125 & full.0$latitude==40.3964996337890625 & full.0$loc_state=="IA"]<-0
full.0$reject[full.0$longitude==-90.41840362548828125 & full.0$latitude==42.261798858642578125 & full.0$loc_state=="IA"]<-0
full.0$reject[full.0$longitude==-88.05780029296875 & full.0$latitude==42.597499847412109375 & full.0$loc_state=="IL"]<-0
full.0$reject[full.0$longitude==-70.8578033447265625 & full.0$latitude==45.26309967041015625 & full.0$loc_state=="MA"]<-0
#full.0$reject[full.0$longitude==-99.8238983154296875 & full.0$latitude==39.378101348876953125 & full.0$loc_state=="KS"]<-0
full.0$reject[full.0$longitude==-71.8114013671875 & full.0$latitude== 43.226398468017578125 & full.0$loc_state=="ME"]<-0
full.0$reject[full.0$longitude==-84.37290191650390625 & full.0$latitude== 46.516201019287109375 & full.0$loc_state=="MI"]<-0
full.0$reject[full.0$longitude==-88.12580108642578125 & full.0$latitude== 45.8069000244140625 & full.0$loc_state=="MI"]<-0
full.0$reject[full.0$longitude==-95.38449859619140625 & full.0$latitude== 49.893100738525390625 & full.0$loc_state=="MN"]<-0
#full.0$reject[full.0$longitude==-87.659698486328125 & full.0$latitude== 39.13330078125 & full.0$loc_state=="IL"]<-0
full.0$reject[full.0$longitude==-92.6667022705078125 & full.0$latitude==43.5 & full.0$loc_state=="MN"]<-0
full.0$reject[full.0$longitude==-95.15219879150390625 & full.0$latitude== 41.146099090576171875 & full.0$loc_state=="MN"]<-0
full.0$reject[full.0$longitude==-96.4860992431640625 & full.0$latitude== 44.339199066162109375 & full.0$loc_state=="MN"]<-0
full.0$reject[full.0$longitude==-96.856903076171875 & full.0$latitude== 46.94329833984375 & full.0$loc_state=="MN"]<-0
#full.0$reject[full.0$longitude==-122.7230987548828125 & full.0$latitude== 45.65000152587890625 & full.0$loc_state=="WA"]<-0
full.0$reject[full.0$longitude==-90.17440032958984375 & full.0$latitude== 38.623699188232421875 & full.0$loc_state=="MO"]<-0
full.0$reject[full.0$longitude==-89.0569000244140625 & full.0$latitude== 36.1074981689453125 & full.0$loc_state=="MS"]<-0
#full.0$reject[full.0$longitude==-86.1322021484375 & full.0$latitude== 41.66500091552734375 & full.0$loc_state=="IN"]<-0
full.0$reject[full.0$longitude==-91.09169769287109375 & full.0$latitude== 32.38970184326171875 & full.0$loc_state=="MS"]<-0
full.0$reject[full.0$longitude==-70.13919830322265625 & full.0$latitude== 43.38330078125 & full.0$loc_state=="NH"]<-0
full.0$reject[full.0$longitude==-70.26080322265625 & full.0$latitude== 43.254199981689453125 & full.0$loc_state=="NH"]<-0
full.0$reject[full.0$longitude==-70.85970306396484375 & full.0$latitude== 43.259700775146484375 & full.0$loc_state=="NH"]<-0
full.0$reject[full.0$longitude==-71.34310150146484375 & full.0$latitude== 42.420299530029296875 & full.0$loc_state=="NH"]<-0
full.0$reject[full.0$longitude==-72.000701904296875 & full.0$latitude== 44.32550048828125 & full.0$loc_state=="NH"]<-0
full.0$reject[full.0$longitude==-114.73889923095703125 & full.0$latitude== 36.0139007568359375 & full.0$loc_state=="NV"]<-0
full.0$reject[full.0$longitude==-114.5 & full.0$latitude==36 & full.0$loc_state=="NV"]<-0
full.0$reject[full.0$longitude==-79.0482025146484375 & full.0$latitude== 43.142398834228515625 & full.0$loc_state=="NY"]<-0
full.0$reject[full.0$longitude==-76.90360260009765625 & full.0$latitude== 43.651401519775390625 & full.0$loc_state=="NY"]<-0
full.0$reject[full.0$longitude==-74.3069000244140625 & full.0$latitude== 45.87689971923828125 & full.0$loc_state=="NY"]<-0
#full.0$reject[full.0$longitude==-88.9839019775390625 & full.0$latitude== 40.53079986572265625 & full.0$loc_state=="IL"]<-0
full.0$reject[full.0$longitude==-82.86150360107421875 & full.0$latitude== 38.648998260498046875 & full.0$loc_state=="OH"]<-0
full.0$reject[full.0$longitude==-120.6968994140625 & full.0$latitude== 45.718898773193359375 & full.0$loc_state=="OR"]<-0
full.0$reject[full.0$longitude==-121.13809967041015625 & full.0$latitude== 45.615001678466796875 & full.0$loc_state=="OR"]<-0
full.0$reject[full.0$longitude==-121.1345977783203125 & full.0$latitude== 45.6139984130859375 & full.0$loc_state=="OR"]<-0
full.0$reject[full.0$longitude==-121.13359832763671875 & full.0$latitude== 45.61370086669921875 & full.0$loc_state=="OR"]<-0
full.0$reject[full.0$longitude==-120.69640350341796875 & full.0$latitude== 45.717998504638671875 & full.0$loc_state=="OR"]<-0
full.0$reject[full.0$longitude==-119.298797607421875 & full.0$latitude== 45.9402008056640625 & full.0$loc_state=="OR"]<-0
full.0$reject[full.0$longitude==-116.9618988037109375 & full.0$latitude== 43.646900177001953125 & full.0$loc_state=="OR"]<-0
full.0$reject[full.0$longitude==-72.59169769287109375 & full.0$latitude== 41.821399688720703125 & full.0$loc_state=="RI"]<-0
full.0$reject[full.0$longitude==-71.5500030517578125 & full.0$latitude== 42.049999237060546875 & full.0$loc_state=="RI"]<-0
full.0$reject[full.0$longitude==-83.0464019775390625 & full.0$latitude== 34.4364013671875 & full.0$loc_state=="SC"]<-0
full.0$reject[full.0$longitude==-96.21970367431640625 & full.0$latitude== 45.221401214599609375 & full.0$loc_state=="SD"]<-0
full.0$reject[full.0$longitude==-96.21779632568359375 & full.0$latitude== 45.221698760986328125 & full.0$loc_state=="SD"]<-0
full.0$reject[full.0$longitude==-84.59220123291015625 & full.0$latitude== 35.440799713134765625 & full.0$loc_state=="SD"]<-0
full.0$reject[full.0$longitude==-84.29560089111328125 & full.0$latitude== 35.166900634765625 & full.0$loc_state=="TN"]<-0
full.0$reject[full.0$longitude==-72.51329803466796875 & full.0$latitude== 42.770599365234375 & full.0$loc_state=="VT"]<-0
full.0$reject[full.0$longitude==-72.44640350341796875 & full.0$latitude== 43.137500762939453125 & full.0$loc_state=="VT"]<-0
full.0$reject[full.0$longitude==-72.3031005859375 & full.0$latitude== 43.6678009033203125 & full.0$loc_state=="VT"]<-0
full.0$reject[full.0$longitude==-70.7082977294921875 & full.0$latitude== 44.311901092529296875 & full.0$loc_state=="VT"]<-0
full.0$reject[full.0$longitude==-88.06999969482421875 & full.0$latitude== 42.016101837158203125 & full.0$loc_state=="WI"]<-0
full.1<-subset(full.0, reject==1, select=-reject)

#Cut controls observations that are also major air polluters
full.1<-full.1[order(full.1$loc_state,full.1$reg_id,-full.1$air),]
full.1$duplicate<-0
full.length<-dim(full.1)[1]
for (k in 2:full.length){
	if (full.1$reg_id[k]==full.1$reg_id[k-1]){full.1$duplicate[k]<-1}
}
full<-subset(full.1,duplicate==0)

#sort the data
#longitude is sorted in reverse; this cuts-down on problems for eliminating later bad segments
full<-full[order(full$loc_state,-full$longitude,full$air,full$latitude),]

############STEP 3: KRIGING WIND DIRECTION###################
#UNIX set directory: cd /Volumes/MONOGAN/
#UNIX installation: R CMD INSTALL CircSpatial_1.0-1.tar.gz
#source: http://cran.r-project.org/src/contrib/Archive/CircSpatial/
#DON'T THINK THIS WORKED IN R: install.packages("/Users/jamie/Downloads/CircSpatial_1.0-1.tar.gz",repos=NULL, type="source")
#source('/Users/jamie/Downloads/CircSpatial 2/R/KrigCRF.R') ####Can I reintegrate this into the source library?

#load revised wind data
wind.0<-read.table('windAngle.txt',header=TRUE,sep="\t")
#View(wind)

#Code radians so that "East" is 0. Essential for the sine and cosine functions.
wind.0$degreesA<- 450-wind.0$degrees
wind.0$oldRadians<-wind.0$radians
wind.0$radians<-pi*wind.0$degreesA/180

#plot the raw data
#png('wind3096.png',width=959,height=593)
us.states<-map("state",fill=FALSE, plot=TRUE, resolution=0)
#points(y=wind$latitude, x=wind$longitude,col='blue',pch=20)
arrow.end.lat<-wind.0$latitude+sin(wind.0$radians+pi)
arrow.end.long<-wind.0$longitude+cos(wind.0$radians+pi)
arrows(x0=wind.0$longitude, x1=arrow.end.long, y0=wind.0$latitude, y1=arrow.end.lat,length=.025,col='blue')
#arrow.plot(wind.0$longitude, wind.0$latitude, u=cos(wind.0$radians+pi), v=sin(wind.0$radians+pi), arrow.ex=0.05, xpd=FALSE, true.angle=TRUE, length=.05, col=1)
#dev.off()

#Create eastings and northings
coordinates(wind.0)<- ~longitude+latitude
proj4string(wind.0)<-CRS("+proj=longlat +datum=NAD83")
#wind<-spTransform(wind.0,CRS("+init=esri:102003"))
#source: http://spatialreference.org/ref/esri/102003/proj4/
wind<-spTransform(wind.0,CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))

#simplify to vectors
x1<-coordinates(wind)[,1]
y1<-coordinates(wind)[,2]
model.direction1<-wind$radians+pi

## Fit An Appropriate Model
FitHoriz1 <- lm(cos(model.direction1) ~ (x1 + y1 + I(x1*y1) + I(x1^2) + I(y1^2))) #2nd order trend
FitVert1 <- lm(sin(model.direction1)  ~ (x1 + y1 + I(x1*y1) + I(x1^2) + I(y1^2))) #2nd order trend
fitted.direction1 <- atan2(FitVert1$fitted.values, FitHoriz1$fitted.values)

## Plot Fitted Model
plot(x1, y1, type="n", asp=1, xlab="", ylab="")
arrow.plot(x1,y1,u=cos(wind$radians+pi),v=sin(wind$radians+pi),col='blue',  arrow.ex=0.05, xpd=TRUE, true.angle=TRUE, length=.05)
arrow.plot(x1, y1, u=cos(fitted.direction1), v=sin(fitted.direction1),  arrow.ex=0.05, xpd=TRUE, true.angle=TRUE, length=.05)

## Compute Residuals
resids1 <- CircResidual(X=x1, Y=y1, Raw=model.direction1,Trend=fitted.direction1, Plot=F)

## Cosineogram 
cosineogram.out<-CosinePlots(x=resids1$x, y=resids1$y, directions=resids1$direction,Lag.n.Adj=1, BinWAdj=1, Plot=TRUE, Cloud=FALSE, Model=FALSE)
abline(h=.24,col=2); abline(v=1000000,col=2)

## Cosineocloud
#CosinePlots(x=resids1$x, y=resids1$y, directions=resids1$direction, Lag.n.Adj=1, BinWAdj=1, Plot=TRUE, Cloud=TRUE)

## Fit of exponential with range=1000000 and sill=.24
#pdf('cosinePlots.pdf')
CosinePlots(x=resids1$x, y=resids1$y, directions=resids1$direction,Lag.n.Adj=1, BinWAdj=1, Plot=TRUE, Cloud=FALSE, Model=TRUE, nugget=0.01,Range=1000000,  sill=0.24, x.legend=.5, y.legend=0.8)
#dev.off()

############STEP 4: INTERPOLATING WIND DIRECTION FOR POLLUTERS###################
air.0<-full

#Create eastings and northings
coordinates(air.0)<- ~longitude+latitude
proj4string(air.0)<-CRS("+proj=longlat +datum=NAD83")
#air<-spTransform(air.0,CRS("+init=esri:102003"))
air<-spTransform(air.0,CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))

#simplify to vectors
x2<-coordinates(air)[,1]
y2<-coordinates(air)[,2]

## Krig to residuals using cosine model. Using previous range and sill estimate.
krig2 <- altKrigCRF(krig.x=x2, krig.y=y2, resid.x=resids1$x, resid.y=resids1$y, resid.direction=resids1$direction, Model = "exponential", Nugget=0.01, Range=1000000, sill=0.24, Smooth=FALSE, Plot=FALSE)
#krig2

## Interpolate Fitted Model
range(x1); range(x2)
range(y1); range(y2)
horiz2<-predict.lm(FitHoriz1,newdata=data.frame(x1=x2,y1=y2))
vert2<-predict.lm(FitVert1,newdata=data.frame(x1=x2,y1=y2))
fitted2<-atan2(vert2,horiz2) #forecasted values

## Plot Estimate Of Direction And Overplot Sample.
points<-length(x2)
map.set<-sample(1:points,1000)
#map.set<-1:points
estimate2=fitted2 + krig2$direction
#pdf('arrowPlot.pdf')
plot(x2, y2, type="n", xlab="", ylab="", asp=1)
arrow.plot(x2[map.set], y2[map.set], u=cos(estimate2[map.set]), v= sin(estimate2[map.set]),arrow.ex=0.05, xpd=FALSE, true.angle=TRUE, length=.05, col="tan")
arrow.plot(x1, y1, u=cos(model.direction1), v=sin(model.direction1), arrow.ex=0.05, xpd=FALSE, true.angle=TRUE, length=.05, col=1)
#arrow.plot(x1, y1, u=cos(wind$radians+pi), v=sin(wind$radians+pi), arrow.ex=0.05, xpd=FALSE, true.angle=TRUE, length=.05, col=1)
#arrow.plot(wind.0$longitude, wind.0$latitude, u=cos(wind.0$radians+pi), v=sin(wind.0$radians+pi), arrow.ex=0.05, xpd=FALSE, true.angle=TRUE, length=.05, col=1)
#dev.off()

#write out data
eastings<-x2
northings<-y2
longitude<-full$longitude
latitude<-full$latitude
angle.full<-estimate2
angle.mod<-fitted2
angle.error<-krig2$direction
air.direction<-cbind(air@data,longitude,latitude,eastings,northings,angle.full,angle.mod,angle.error)

write.csv(air.direction, "anyPowerPlantDirection.csv",row.names=FALSE)

############STEP 5: DIRECTION TO LEEWARD BORDER###################
#relabel
#full<-read.csv('anyPowerPlantDirection.csv')
full<-air.direction

#load map polygon data
map.states<-readShapePoly(fn="tl_2011_us_state.shp")

#set up input and output objects
states<-c('AL','AR','AZ','CA','CO','CT','DE','FL','GA','IA','ID','IL','IN','KS','KY','LA','MA','MD','ME','MI','MN','MO','MS','MT','NC','ND','NE','NH','NJ','NM','NV','NY','OH','OK','OR','PA','RI','SC','SD','TN','TX','UT','VA','VT','WA','WI','WV','WY')
numbers<-c(41,53,6,5,1,15,11,44,46,49,14,16,29,24,25,52,51,47,45,28,8,54,9,36,22,40,3,43,56,33,7,13,21,30,27,34,32,35,48,42,2,50,20,12,37,38,55,4)
factories<-rep(NA,48)
in.polygon<-rep(NA,48)
no.polygon<-rep(NA,48)
no.segments<-rep(NA,48)
deviation.long<-rep(NA,48)
deviation.lat<-rep(NA,48)
distances<-c(0)

#loop to create ppp objects and create distance measures
for(i in 1:48){
	#i<-1 #index for checking loop code one item at a time
	assign(states[i],SpatialPolygons(Srl=list(map.states@polygons[[numbers[i]]])))
	data.for.state<-subset(full,loc_state==states[i])
	factories[i]<-dim(data.for.state)[1]
	
	#necessary information for window creation
	no.polygon[i]<-length(get(states[i])@polygons[[1]]@Polygons)
	full.coords<-list()
	loc.min.x<-rep(NA,no.polygon[i]); loc.max.x<-rep(NA,no.polygon[i]); loc.min.y<-rep(NA,no.polygon[i]); loc.max.y<-rep(NA,no.polygon[i])
		
		#create window features	
		if (no.polygon[i]==1) {
		state.trim<-as.matrix(get(states[i])@polygons[[1]]@Polygons[[1]]@coords[-1,])
		state.bound<-owin(xrange=c(min(state.trim[,1]),max(state.trim[,1])),yrange=c(min(state.trim[,2]),max(state.trim[,2])),poly=list(x=rev(state.trim[,1]),y=rev(state.trim[,2])))
		} else {
			for(j in 1:no.polygon[i]){
				state.trim<-as.matrix(get(states[i])@polygons[[1]]@Polygons[[j]]@coords[-1,])
				loc.min.x[j]<-min(state.trim[,1]); loc.max.x[j]<-max(state.trim[,1]); loc.min.y[j]<-min(state.trim[,2]); loc.max.y[j]<-max(state.trim[,2])
				full.coords[[j]]<-list(x=rev(state.trim[,1]),y=rev(state.trim[,2]))
				}
		min.x<-min(loc.min.x); max.x<-max(loc.max.x); min.y<-min(loc.min.y); max.y<-max(loc.max.y)
		state.bound<-owin(xrange=c(min.x,max.x),yrange=c(min.y,max.y),poly=full.coords)
				}

	#combine window feature into ppp object
	state.ppp<-ppp(x=data.for.state$longitude, y=data.for.state$latitude, window=state.bound, marks=as.factor(data.for.state$air))
	assign(paste(states[i],'ppp',sep='.'), state.ppp)
	in.polygon[i]<-state.ppp$n
	#attr(state.ppp,'rejects')$x
	#attr(state.ppp,'rejects')$y
	#plot(state.ppp)
	#full[full$loc_state=='AL',c('latitude','longitude')]
	
	#generate lines for distance measure, create psp object
	reach.x<-state.ppp$x+10000*cos(data.for.state$angle.full)
	reach.y<-state.ppp$y+10000*sin(data.for.state$angle.full)
	#plot(KY.ppp)
	#arrows(x0=state.ppp$x,x1=reach.x,y0=state.ppp$y,y1=reach.y,length=.05,col='blue')
	initial<-psp(x0=state.ppp$x, y0=state.ppp$y, x1=reach.x, y1=reach.y, window=owin(xrange=state.ppp$window$xrange, yrange=state.ppp$window$yrange),check=FALSE)
		#window=owin(xrange=c(-10000,10000),yrange=c(-10000,10000)))#,
	assign(paste(states[i],'psp',sep='.'), initial[,state.bound])
	no.segments[i]<-dim(get(paste(states[i],'psp',sep='.'))$ends)[1]
	angles<-angles.psp(get(paste(states[i],'psp',sep='.')))
	
	#drop duplicate lines from the PSP for the same point
	#count how many times a segment overlaps with the state border
	keep<-rep(1,no.segments[i])
	splice.count<-rep(NA,length(reach.x))
	for(j in 1:length(reach.x)){
		temp.psp<-psp(x0=state.ppp$x[j], y0=state.ppp$y[j], x1=reach.x[j], y1=reach.y[j], window=owin(xrange=state.ppp$window$xrange, yrange=state.ppp$window$yrange),check=FALSE)
		splice.count[j]<-nsegments(temp.psp[,state.bound])
	}
	mult.segs<-which(splice.count>1)
	overshoot<-c(0)	
	add.overshoot<-	splice.count[which(splice.count>1)]-1
	if(max(splice.count)>1 & length(mult.segs)==1){ 
		starts<-mult.segs+overshoot+1
		stops<-starts+add.overshoot-1
			for(k in 1:length(starts)){keep[(starts[k]):(stops[k])]<-0}
	}
	if(max(splice.count)>1 & length(mult.segs)>1){ 
		for(k in 2:length(mult.segs)){overshoot[k]<-overshoot[k-1]+add.overshoot[k-1]}
		starts<-mult.segs+overshoot+1
		stops<-starts+add.overshoot-1
			for(k in 1:length(starts)){keep[(starts[k]):(stops[k])]<-0}
	}

	#compute the haversine distance
	segs<-get(paste(states[i],'psp',sep='.'))$ends[keep==1,]
	#as.numeric(segs$x0==data.for.state$longitude) #nuisance segment checker
 		#eliminate 4 nuisance segments that don't belong:
 	#	if (i==15){segs<-segs[-488,]}
	# 	if (i==30){segs<-segs[-227,]}
 	#	if (i==36){segs<-segs[-540,]}
 	#	if (i==40){segs<-segs[-396,]}
	distances<-append(distances,haversine(segs$x0,segs$x1,segs$y0,segs$y1))

	#data checking
	no.segments[i]<-dim(segs)[1]
	deviation.long[i]<-sum(data.for.state$longitude-segs$x0)
	deviation.lat[i]<-sum(data.for.state$latitude-segs$y0)
}

###checking  that sample sizes align and that segments properly correspond to points###
#table(factories[]-in.polygon[])
#table(no.segments[]-in.polygon[])
#table(deviation.long); table(deviation.lat)
#length(distances)
#dim(full)

#add downwind distances to main data
full$down<-distances[-1]

#standardized downwind
full$sDown<-full$down/ave(full$down,full$loc_state,FUN=max)

############WRITE OUT THE FILE###################
write.csv(full, 'anyPowerPlantDist.csv', row.names=FALSE)
